library(tidyverse)
library(ggstance)
library(magrittr)
library(readstata13)
library(lmtest)
library(sandwich)

d1 <- "C:/Dropbox/NyhanPorterReiflerWood/replication data and code/Study 1/pooled-data-for-R.dta" %>% 
  read.dta13() %>% 
  tbl_df %>% 
  select(
    crimechange,
    stataccuracy,
    articleaccuracy,
    trumpsupport,
    trumpcrime,
    trumpcrimeXtrumpsupport,
    correction,
    correctionXtrumpsupport,
    questionstats,
    questionstatsXtrumpsupport,
    fbict,
    fbictXtrumpsupport,
    clintonsupport, trumpsupport, turk
  ) %>% 
  filter(
    turk == 0 &
      (
        trumpsupport == 1 |
          clintonsupport == 1
      )
  ) %>% 
  mutate(
    id = 1:n()
  ) %>% 
  # generate summary factors
  mutate(
    cond = case_when(
      trumpcrime == 1 ~ "trumpcrime",
      correction == 1 ~ "correction",
      questionstats == 1 ~ "questionstats",
      fbict == 1 ~ "fbict",
      TRUE ~ "birdwatch"
    ) %>% 
      factor(
        c("birdwatch",
          "trumpcrime",
          "correction",
          "questionstats",
          "fbict")
      ),
    vote16 = case_when(
      trumpsupport == 1 ~ "trump voter",
      clintonsupport == 1 ~ "clinton voter"
    ) %>% 
      factor(
        c("clinton voter",
          "trump voter")
      )
  ) 

lm1 <- lm(
  crimechange ~ cond * vote16, 
  data = d1
)

# we recover the Table 1 estimaes

lm1 %>%
  coeftest(
    vcov = vcovHC(lm1, type="HC1")
  ) %>% 
  round(2)

eg1 <- expand.grid(
  cond = d1$cond %>% levels %>% rev,
  vote16 = d1$vote16 %>% levels
)

eg1 %>% 
  mutate(
    fit = lm1 %>% 
      predict(
        eg1
      ),
    se = lm1 %>% 
      predict(
        eg1, se.fit = T
      ) %>% 
      use_series(se.fit),
    lo = fit %>% 
      subtract(
        se %>%
          multiply_by(1.96)
      ),
    hi = fit %>% 
      add(
        se %>%
          multiply_by(1.96)
      ),
    lab = fit %>% 
      round(2) %>% 
      str_pad(width = 4, side = "right", pad = "0")
  ) %>% 
  filter(
    cond %>% 
      equals("birdwatch") %>% 
      not
  ) %>% 
  ggplot() +
  geom_pointrangeh(
    aes(
      xmin = lo, xmax = hi, x = fit, y = cond, linetype = vote16
    ),
    shape = 21,
    fatten = 20,
    fill = "grey90"
  )  +
  geom_text(
    aes(fit, cond, label = lab),
    size = 3
  )

lm2 <- lm(
  stataccuracy ~ cond * vote16, 
  data = d1
)

lm3 <- lm(
  articleaccuracy ~ cond * vote16, 
  data = d1
)


eg1 %>% 
  mutate(
    fit = lm2 %>% 
      predict(
        eg1
      ),
    se = lm2 %>% 
      predict(
        eg1, se.fit = T
      ) %>% 
      use_series(se.fit),
    lo = fit %>% 
      subtract(
        se %>%
          multiply_by(1.96)
      ),
    hi = fit %>% 
      add(
        se %>%
          multiply_by(1.96)
      ),
    lab = fit %>% 
      round(2) %>% 
      str_pad(width = 4, side = "right", pad = "0"),
    dv = "stat accuracy"
  ) %>% 
  bind_rows(
    eg1 %>% 
      mutate(
        fit = lm3 %>% 
          predict(
            eg1
          ),
        se = lm3 %>% 
          predict(
            eg1, se.fit = T
          ) %>% 
          use_series(se.fit),
        lo = fit %>% 
          subtract(
            se %>%
              multiply_by(1.96)
          ),
        hi = fit %>% 
          add(
            se %>%
              multiply_by(1.96)
          ),
        lab = fit %>% 
          round(2) %>% 
          str_pad(width = 4, side = "right", pad = "0"),
        dv = "article accuracy"
      )
  ) %>% 
  filter(
    cond %>% 
      equals("birdwatch") %>% 
      not
  ) %>% 
  ggplot() +
  geom_pointrangeh(
    aes(
      xmin = lo, xmax = hi, x = fit, y = cond, linetype = vote16
    ),
    shape = 21,
    fatten = 20,
    fill = "grey90"
  )  +
  geom_text(
    aes(fit, cond, label = lab),
    size = 3
  ) +
  facet_grid(dv ~ .)


